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Abstract. We show how to exploit algebraic relations among the operators (or 
matrices) which constitute the non-equilibrium matrix product steady state of a 
boundary driven quantum spin chain to find partial differential equations determining 
all the m-point correlation functions in the continuum (or thermodynamic) limit. 
These partial differential equations, the order of which is determined by scaling of the 
non-equilibrium partition function, are readily solved if we also know the boundary 
conditions. In this way we can avoid resorting to representation theory of the matrix 
product algebra. We apply our methods to study the distributions, or moments, of 
the magnetization and the spin current observables in boundary driven open XXZ spin 
chains, as well as some connected correlation functions. Interestingly, we find that the 
transverse connected correlation functions are thermodynamically non-decaying and 
long-range at the isotropic point A = 1. 
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1 . Introduction 

Matrix product ansatz (MPA) for steady states has a long history in non-equilibrium 
physics. MPA was hrst used to hnd the non-equilibrium steady states of a ID 
asymmetric exclusion model analytically [T] and were also extended to other classical 
driven diffusive systems [2]. Later MPA has been also applied to steady states of open 
quantum spin systems, which are described as fixed points of Lindblad master equations 
[3]. In these models the spin systems are coupled to two Markovian baths, which drive 
the system out-of-equlibrium. These types of setups have attracted lots of attention 
lately in the context of transport theory (see, e.g. Refs. [U El El 0 El El Eli E2l ESI 
ESIES]), as their study has become accessible to experiments through various quantum 
simulations techniques, such as those using cold atoms EEl EZl EHl ES]. They also find 
potential application in the context of e.g., quantum computing [201 EH 122]. 

When we are studying the statistical properties of such systems either the 
continuum or thermodynamic limit is the most interesting. Due to the richness of 
phenomena and highly non-trivial nature of out-of-equilibrium systems interesting and 
useful information can be gained from studying not merely the expectation value of 
physical quantities, but their statistical properties (fluctuations) as well. For instance, 
one may be interested in e.g., cumulants, correlators, connected correlation functions in 
various contexts. 

Calculation of the statistical properties can be greatly eased through the use of large 
deviation theory and the related method of full-counting statistics, by means of which 
one can in principle calculate the full probability distributions of physical quantities 
of interest [23] • These two methods were only very recently applied in full to quantum 
systems m and even more recently to both non-interacting (see e.g., [251 EH] EH EH] ) ^md 
interacting [291 EDI [31] many body systems. Using these methods open quantum systems 
can be also seen to exhibit interesting properties near phase transitions [32l l33l EH EH] 
and one may also extend them to study closed systems [361 EZ] • 

However there are only a few known exact analytical results for the full counting 
statistics of non-interacting many-body quantum systems [261 EH EH, and likewise, 
there is only a single analytical result for an interacting case of the XXZ spin chain 
[29] . The latter was achieved only perturbatively in system-bath coupling. In fact, 
analytically studying interacting many-body quantum systems under the open quantum 
system framework seemed like a formidable task, but became feasible after two recent 
results [Hi for the open XXZ spin chain, later understood through the underlying 
quantum integrability of the system [TUI EHl EH] (for a review see [I])- 

Throughout our work we shall call the matrices constituting the matrix product 
steady state as auxiliary space operators (ASO). In the aforementioned solutions the 
ASOs fulhl certain algebraic relations, which we will call the matrix product algebra. 
To compute observables one must usually employ an appropriate representation of the 
algebra. We will show however that what one requires, in principle, are only the dehning 
relations of the algebra, which in the continuum limit (or thermodynamic limit) lead 
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to partial differential equations (and of course the corresponding boundary conditions 
we need to solve these equations). What we study in this article is essentially a simple 
generalization of the procedure employed in [8] to calculate the magnetization profiles 
and 2-point connected longitudinal spin correlation function (the same method was also 
later used in Ref. [ID] to compute the profiles and currents - but not the correlations - 
for more general, so-called twisted boundary conditions). 

Using this method we then find explicit expressions for several m-point connected 
correlation functions and spin current fluctuations for the boundary driven open XXZ 
spin chain Interestingly, for the critical A = 1 case the transverse connected 

correlation functions are non-decaying and thus exhibit genuine long-range order, similar 
to what has been previously observed numerically in a related case m We also study 
the probability distribution of total magnetization and the moments of the spin current 
operator. 

Note that we compute these correlators for the non-equilibrium steady state, which 
can be contrasted with other dynamical studies of both open quantum systems w 
and systems undergoing quantum quenches |13]. We should also note that Verstraete 
and Cirac [H] introduced continuous matrix product states (cMPS) for quantum fields. 
Our approach is not related to this. Instead of constructing matrix product states for 
quantum fields, which are continuum limits of lattice theories, we will take a discrete 
matrix product and study the continuum limit. 

In this paper we discuss a general procedure for computing the continuum limit of 
a steady state (assumed to be given in the form of a matrix product state). A key step 
in taking this continuum limit is a perturbative expansion in lattice spacing, the validity 
of which is not known. The second result, is the computation of connected correlators 
and fluctuations of current in the steady state for the open maximally driven XXZ spin 
chain in the continuum limit. Using a known discrete solution for this model [H] we can 
check the validity of our method. 

More specifically, in Sec.j^we review the properties of matrix product steady states. 
In Sec. [^we outline our method for computing the continuum limit of the steady state 
equation and steady state (under certain assumption discussed). Later, in Sec. [fusing 
this method (and also aided with the known solution for the discrete steady state [H]) 
we compute the correlation and connected correlation functions for the steady state of 
the aforementioned open maximally driven XXZ spin chain, focusing mostly on the non¬ 
trivial isotropic XXX case. Afterwards, in Secs. andwe study the fluctuations of spin 
current and total magnetization in the steady state of the open maximally driven XXZ 
spin chain, aided by our previous computation of the connected correlation functions. 

2. Matrix product steady states 

We will be interested in non-equilibrium steady states (NESS) poo of one-dimensional 
spin-1/2 systems (quantum spin chains). Let the system have n sites described by 
a 2”-dimensional Hilbert space ^ on which act operators constructed from the Pauli 


Connected correlations, fluctuations and current of magnetization in the steady state.. A 


matrices, af, cr|, <7° := 1, where j = 1,... ,n labels the site position. Let the dynamics of 
the system, described by the density matrix p{t), be determined by a qnantum Lionville 
eqnation, 

^p(() = (1) 

where can be nnderstood as a snperoperator acting on the space of operators 
e^(^), spanned by the Panli matrices. The space may also be considered 

as a Hilbert space itself if one defines an inner prodnct in the Hilbert-Schmidt sense, 
i.e., A,B e ^(^), {A,B) = iiA^B. From Eq. 0 the defining eqnation for the NESS 
is, 

^Poo = 0 , ( 2 ) 


The key assnmption we will nse is that the NESS is given in the form of a homogenous 
MPA, 


poo 


tlSn 


(3) 


where. 


= (L|L®"|R), 


(4) 


such that. 


Sn = {L 


/o. 

0 -] 

' 0, 

02 / 


|R), 


(5) 


where Oj G End(6s) are called the auxiliary space operators (ASO)[|] acting over the 
vector space G^, which is also the space on which the representation of the symmetry 
algebra of the model acts. Importantly, these operators satisfy some algebraic relations, 
which are assumed to be known. The states |R), |L) G G^ are referred to as the boundary 
vectors. We will also define four important operators, Oq, O^, Ox, and Oy, 


Oq Oi + O2, Oz Oi — O2, (6) 

0x:=0+ + 0_ Oy ;=-i(0+-0_) (7) 


One usually assumes that the representation (used to calculate the steady state 
explicitly) of the algebra satisfied by the auxiliary space operators is known. We will 
not do so here, but will first illustrate the approach one takes if it is known. Observables 
can be calculated from the MPA in the following simple way, provided one knows the 
representation of the ASO algebra. 

Define a general, not necessarily local, operator = cr“W 2 ^... a"", where 

1.. .n label the sites and aj G {x,y,z,0} denote the components of the corresponding 
Pauli matrices use]. Its expectation value in the steady state is given by (from Eq. 
0 ). 


f These operators were sometimes also called vertex operators in Refs. mzii, but they have no 
relation to the more standard concept of vertex operators (see e.g., [75] 1. 
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An object which will be of central importance is the so-called non-equilibrium partition 
function, 

Xi = trS'n (9) 

It is related to currents flowing through the system in the NESS for a wide variety of 
exactly solvable ID systems, including both classical processes [2] and various out-of- 
equilibrium quantum spin chains [1]. 

Using the MPA form in Eq. ([^ Eq. ([^ can be written as. 


(L| 


{B, 




.) = 


tr. 



ai,...an 


|R) 


trS'n 


( 10 ) 


where we have used the fact that the trace tr := trp is taken only over the physical 
2”-dimensional Hilbert space and not over the auxiliary space by definition (Eq. 
(|^). Then it is merely a matter of simple matrix multiplication (in the physical space) 
and using repeatedly the property of the trace that tr(A 0 5) = tr(A)tr(5) (together 
with definitions Eqs. ([^, ([^) to find that, 

, (L|0^,0^,...0^JR) 

^ (L|0^|R) ’ ^ ^ 


where aj G {x, y, z,0}. Note that this also gives us that the non-equilibrium partition 
function (j^ can be written as. 


= tr^„ = (L|0((|R). 


( 12 ) 


We can thus define a mapping from expectation values of observables to their 
corresponding auxiliary space operators, e.g.. 


{(^> 1 ) = 


(L|0' 




|R) 


(L|OS|R) 


(13) 


In order to actually calculate the expectation values of operators in the NESS 
we also need know the representation of the ASOs Oj. For interacting systems these 
representation are generically inhnite dimensional. Even though the representations are 
near-diagonal in the integrable (in the sense of Ref. i) cases and thus allow for efficient 
computation, calculating the expectation values when the operators are not ultra-local 
(acting only on a single site) or for very large systems is impossible. Using the approach 
discussed in the next section we show that we can bypass this difficulty (at least up 
to a multiplicative prefactor) by employing only the asymptotic (n —)■ cx)) form of the 
non-equilibrium partition function together with the algebra satisfied by the ASOs, 
to calculate all the m-point correlators (and equivalently the entire steady state of the 
system) in the continuum limit without resorting to an explicit representation of the 
ASOs. 
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3. Continuum limit of the NESS 

Our procedure is similar to that used for one and two-point functions in [S] and later in 
Motivated by these known examples we will consider only ASOs which satisfy at 
most cubic algebraic relations of the form, 


+ fcg 2O0O0O0 -|- ^3 3O0O0OQ 

-|-/C2,lOoOo + ^2,20 oOq, -|- = 0, (14) 


where a G {x, y, z} for some constants kfy In principle our approach can be generalized 
to other cases as well, though we will not discuss this here. Let us now introduce a 
lattice spacing a = 1/n, such that the total length of the system is unity, and so the 
continuum limit a —)■ 0 corresponds to the thermodynamic limit n —)■ cx). 

We wish to hnd a set of partial differential equations for a m-point correlator, 

"), where aj G {x,y,z} for each of the m operator 




laflaT 

' Jl 32 


■ 

Jn 


coordinates. 

We will first find them for ai (ji) by multiplying Eq. 


'“'o '^a 2'^0 




0„3 . . . O, 


3m 


(14) for a = «!, by 
|R) from the right and by (L|0;^^~^ from the left 
and divide it by the non-equilibrium partition function ^ = (L|Og|R). We then use 
Eq. ( [II| to find similarly to jlUl 0] , 


7 Qi i iOL\ , 7 Q;i . ,am 

^3,2^ji+l,j2+3,...j^+3 ^3,3^j 


■^jl+ 2 ,J' 2 + 3 ,.--Jm +3 


I I 7 ai ^a:i,...,Q;Tn;n. \ ^n-1 ^n-2 


‘=^'n 


= 0, (15) 


where the superscript n over denotes that the this correlator is computed for 

system size n. Define Xk = Xk + a ■.= {jk + ^)/'n, with inter-site (lattice) spacing 
denoted as a := 1/n, and likewise C°‘^’-°‘"'{xi... Xm) ■= Note that like in 

[H m 00] one may instead equivalently define Xk = [jk ~ l)/(’^ ~ 1) and a = l/(n — 1) 
as we will do in the next section. Assume that we can expand for large n as. 


6^ 

‘=^n—l 

6^ 


oo 



m=0 


(16) 


where first few coefficients z^^\ etc, may be vanishing. As we will see later, this type 
of expansion can be performed for the open maximally driven XXZ spin chain at A < 1 
and for some other models such as open S't/(A)-symmetric quantum gases [39|, XXZ 
spin chains with twisted boundary driving (also for A < 1) [TUI 00], the open spin-1 
Lai-Sutherland chain [38] and some other cases discussed in the review article [1]. We 
then take the continuum limit n —)■ cx) by expanding in l/n[§| 

= ++ = (17) 

n n 


§ Note that we have suppressed the superscript ai... in ... Xm)', when we do this we refer 

to a general m-point correlator, i.e., C(xi,... x^) = (aii,... x^). 
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and performing Taylor series expansioE[j|] in a = 1/n around 0 to find 
system of partial differential equations up to arbitrary order in 1/n (where, 
kfl ■= kff + + kfl, kfl := kf\ + and ...,Xm):= C^’^\x)), 


'^3,2 


''2,1 


'' 2,2 


an infinite 
for brevity, 


+ kf^)C^^\x) =0, 

+ iJ')C''>(2) + (k-l, - k« + \{kf )xi - i“;])4.C<»>(2))+ 

- 1)4,C'<'»(2) = 0, 

(^( 2 )^( 0 ) (-) ^ ^( 1 )^( 1 ) (^)) + 

- I)"j''>ia2,c>''>(2)} + (k-gx, -1) + i«:e0z'">s..c«>|(2)+ 

(*^3,3 - ^“l + - *:J,i)3"”)8„C<"(2) + 

^(^ 2,1 ( 2^1 - 1 ) + kf^2Xi){xj - l)z^°'>d,,,d,,.C^^\x)+ 

^Ki + + (*'a(l - 2x0 + krx^,)z«»]alc«»(x) = 0 

: (18) 


We then continue by deriving the equations for j 2 by hrst multiplying Eq. 0 by 

from the left, etc. and likewise for all aj. We are hnally left with 
a set of coupled partial differential equations for every order in 1/n. These are in general 
complicated for arbitrary orders, but if we only focus on the leading order (l/n)° they 
are generically quite simple. The leading order will be determined by the hrst non-zero 
equation in the set Eq. (18). For instance if z^^'^kfl + kf^ flz 0 the leading order is given 
by C(°)(f) = 0. 

Note that everything, except the boundary conditions, is fully determined by the 
algebraic relations Eq. (14) and the asymptotic scaling n —)■ cx) of the non-equilibrium 
partition function ^ Eq. (§. 

We did not consider the representation of the ASOs at all. One may object that 
the representation is relevant when one wishes to hnd the boundary conditions to solve 
these partial differential equations and that it comes into play via the boundary vectors 
|R) and |L), which we used when deriving Eq. (15). However, in the leading order at 
least this can be circumvented for a quite general set of Liouvillians Eq. ([^ which dehne 


II We assume that we can perform this expansion. The points when the indices coincide can introduce 
extra boundary conditions and even cause our expansion to fail when these indices are close to each 
other in the continuum limit. For the open maximally driven XXZ case we study later we have the 
added benefit of having a discrete solution for the NESS which was previously found to check our results 

0. 
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the non-equilibrium steady state Eq. (|^, 

Poo = (-^0 + -=^1 + ■^n)Poo = 0 (19) 


where acts in general in the bulk of the system and and act only ultralocally 
(on the boundary sites 1 and n, respectively). These are the types of Liouvillians one 
most often encounters when studying non-equlibrium matrix product steady states. 

As mentioned previously, we work with one-dimensional spin-1/2 systems and thus 
we take that is associated with a model with hnite lattice spacing a. For instance, 
in the case of the Lindblad master equation (which we will study in more detail in the 
next section) where H is the Hamiltonian of the system. 

Assume that can be written in terms of local two-site interaction operators, 
where j, j + 1 denotes the sites on which the operator acts. We then 
perform the continuum limit as before by hrst setting a = 1/n, x = j/n, x+a = {j+l)/n. 
Namely, 


^ i{x = ^) 1 ^ i{x + a = 


( 20 ) 


Formally then, when taking the continuum limit, a = l/n—)-0we may expand for small 
a, =^o(a) = -^0°^ + ^(a), where the a in J:fo(a) denotes that we are now dealing with an 
operator which depends on lattice spacing a after we used Eq. (20). 

Likewise we may formally take the continuum limit of the NESS by hrst writing 
out in the operator basis. 


Poo = at I 1 ^ ^ 

k^m 


( 21 ) 


k a 


where is a normalization coefficient (such that trpoo=l), and then taking the same 
continuum limit Poo(a) = P^^^ + p*'^^a+ ^(a^). We do this in the following manner. First 


rewrite Eq. (21) as discussed. 


poo 


\ k Oi 


{x = ka))a°'{x = ka) 


+ (o'“((Ei = ka)a^{x2 = ma))a°'{xi = ka)a^{x2 = ma) ... j . ( 22 ) 

k^m a,0 


We then formally expand the correlation functions in Eq. (22) for small a = 
1 /n, {a°‘{x = ka) = {a‘^{x))^^'’ + {a^{x))^^'^a + {a^{xi = ka)a'^{x2 = ma)) = 

{a°'{xi)a^{x2))^^^ + ((T"(a;i)(T^(a;2))*'^^a + C’{a‘^),... Plugging this back into Eq. (22) 
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and gathering terms in a we identify, 


p(o) = 




(x = fca))*'°V"(x = ka) 


+ lim ((x“(xi = ka)a^{x 2 = ma))^°V“(xi = ka)a^{x 2 = ma)... 


a =-^^0 


(23) 


p<'> = N 


k^m q ,/3 

lim 

a=i-s>0 

k n 


= A;a))*'^^cr"(x = fca) 


v(l) 


lim y y {a°‘{xi = ka)fl\x 2 = ma)) V“(xi = ka)m{x 2 = ma)... 

a=--5'0 , , „ 

k^m a,13 


(24) 


Let ns pause to make a few comments. When the difference between xi and X 2 , etc. is 
of the order of lattice spacing a the expansion may be ill-dehned. In fact, this will turn 
out to be the case for the A < 1 case latter studied. When we know the discrete solution 
for the NESS, which will be the case when we will later study the open XXZ spin chain, 
we can use this solution to check our results for the continuum limit. Otherwise, one 
may have to simply assume that the continuum limit can be taken. Furthermore, the 
normalization coefficient N depends on a = 1/n. However, it can be seen to cancel in 


the equation for the NESS (19) and thus does not influence the physical results. 


We also assume that Afi and do not depend on a. We then have in the leading 
order a°, 

+ X + (25) 

where is essentially almost equivalent to knowledge of all the correlators in the 
continuum limit. 

Note that this continuum limit is the same as the one we took when calculating the 


differential equations for the correlators Eq. (18). In other words, if we already know 


the discrete solution for the NESS poo taking the continuum limit for the correlators as 
we did when Ending Eq. (18) also gives the solution perturbatively in lattice spacing 


a. Since Afi and are ultralocal (assumed to be acting on one site each) solving Eq. 


(25) and thus obtaining the boundary conditions needed to solve the leading order of 


the set of partial difierential equations in Eq. (18) is simple. 


It is important to note that even though we may circumvent the issue of not knowing 
the boundary vectors |R) and |L) and the representation of the ASOs using the above 
discussed approach, it is not necessary to do so. In case of the already solved problem 
of the steady state of the open non-equilibrium boundary driven XXZ spin chain [8] in 
terms of an MPA the boundary vectors and the representation of the ASOs are known. 
One can then simply use this to find the appropriate boundary conditions when solving 


the partial differential equations (18) in a manner similar to what was done in [SI SO] 
for a less general set of correlators. 


We will use Eq. (25) in the next section to find the boundary conditions for the 


partial differential equations determining the NESS of a boundary driven open XXX 
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spin chain. 

We will now turn to an example of a previously solved MPA steady state of an open 
non-equilibrium boundary driven XXZ spin chain and compute the m-point correlators 
in this case. 

4. The m-point correlators of the maximally boundary driven XXZ spin 
chain 

The Lindblad master equation is a useful tool for describing out-of-equilibrium physics 
m It can represent both driving and decoherence by a set of infinite baths coupled to 
a system under the Born-Markov and rotating wave approximations [3]. It also has an 
important property of being the most general form of a time-local Markovian quantum 
master equation which is both completely positive and trace preserving. The Lindblad 
master equation is, 

= ^p{t) = + 9(p(t)). 

9(p(t)) [upmi - 1 {dn.p(«)}). (26) 

k ^ ' 

where we will take H to be the XXZ spin chain Hamiltonian, 

n 

^ = ^7+1 + ^i+i ’ (27) 

i=i 

and Lindblad operators acting only on the boundary sites 1 and n, 

Li = y/eaf, Ln = \/ea~. (28) 

They represent maximum driving, that is the left bath decoherently only injects 
magnetization into the system and the right one only takes it out. The equation for the 
NESS is exactly solvable (where we define the superoperator (adLr)p = [H,p]), 

■^Poo = -iadiLpoo + ^(Poo) = 0, (29) 

due to the underlying integrability structure of the XXZ spin chain. The solution was 
found in [8] (see jl] for a more comprehensive overview). 

f.l. The isotropic point A = 1 

At the isotropic point A = 1, it is known that the AS Os satisfy the following cubic 
algebraic relations jH HOI 136] . 

[Oo, [Oo, 0„]]-F 2{Oo, 0„} - 8p^O„ = 0, a = x,y,z, ^ (2'^) 

It is also known that for e S> 1/n the partition function scales as [3], 
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Using the method discussed in the previous section we immediately arrive at a set of 
decoupled second order partial differential equations for the m-point correlators in the 
continuum hmil|^ 


dxl 


C^^\x) =k = l,...,m, 


(32) 


where we emphasize again to avoid confusion, that denotes the leading order in 

1/n of a general m-point correlator C{x), 


C{x) := 


cy .2 
31 ^32 


3n 


:= {a'^flx,)a^flx2)...a'^-{x^)), 


i.e, 


C{x) = ^C^’^\x) 


fc =0 


Xk = (33) 

n 


(34) 


Now we will show how one can obtain the boundary conditions as outlined briefly 
in the previous section, Sec. First we take x = j/n = ja and set af —cr“(x), 
—)■ a'^(x + a) in the Hamiltonian (27) (for A = 1), while keeping the length 
hxed na = 1, and then expand in lattice spacing a. To do this hrst look at a pair 
local densities = X)a=x y zdiscrete Hamiltonian (27) (noting that 


H = which in the continuum limit go 


as, 


(35) 


hjj^i + hj-ij h{x,x + a) + h{x — a,x) = cr"(a; —a)(j"(a;)-|-(T“(a;)(j"((r-f a). 

Q:=x,y,z 

Expanding in a, we hnd, 

02 Q/ \ 

h{x,x + a) + h{x — a,x) = ^ 2cr"(a;)cr"(a;)-|-2cr"(x)— ^ -|-^(a^), (36) 


a=x,y,z 


i.e., the term of order a containing the hrst derivative cancels and the hrst non-trivial 
term is the one with the second order derivative. It may be interesting to note that this 
cancellation of the terms with the hrst derivatives mimics the one we got when deriving 
Eq. (32). Now using the (Riemann) sum dehnition of the integral and the basic algebra 
of the Pauli matrices we arrive that in the leading orders of a the Hamiltonian is. 


H = -6 -I- a 

a 


jy-z 


Q(=x,y,z 




(37) 


The leading term here in Eq. (37) is divergent as a —)■ 0. However, this is not much of a 
problem as this term commutes with any operator in ad H in the steady state equation 
(29) and thus has no inhuence on the hnal result for the NESS. Likewise the term of 
order commutes with any operator and does not inhuence ad 77. 


^ Note that it turns out due to the cubic algebra of this problem Eq. (301 that the leading order terms 
after performing the expansion via Eq. (181 do not depend on e and p. 
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Therefore, the hrst non-trivial term in Eq. (37) is of order a. We note for the 
interested reader that, after using integration by parts, this term corresponds to what 
is known as the quantum Landau-Lifshitz model (e.g., |1HI 09]) or the SU{2) quantum 
continuous Heisenberg magnet |50j . 

We can also expand poo using Eq. (21) for small a. Then the leading order in a is 
given by (using Eqs. (25), (23) and Eq. ([2^), 

= 0, (38) 


where ^ is given by Eq. (26) and Li = Ln = (1). Note that assuming 

that £ 7 ^ 0 we can cancel it in Eq. (38). Splitting the dissipator into ^ 
where ^ p| j, we observe the following action on the 

basis operators. 


^i(l) =£a^(0). 


k(cTflO)) = -eaflO), 


^.(cr^’^(l)) = -|u^’^(l). (39) 


Using this Eq. (39) and requiring a solution to Eq. (38) for each of the operators in the 
basis Eq. (21) we arrive to the following boundary conditions, 

C^°n^)Uqxfc=o) = = 0 )cr"'=+i (0:^+1)... cr“™(a;^)) = 

(a^flx,)... ... a^-(xj}, (40) 

C^°n^)Uqxfc=i) = (u"Hxi). .. a^^-flxk-i)o-flxk = l)o-^'°+flxk+i). . .cr“™(a;m)) = 

- (a“i(a;i)... a^'‘-flxk-i)o-'"'^+flxk+i)... a^’^(xm)}, (41) 

= (o-'^flxi)... a°‘^-flxk-i)o-'^’^(xk = 0, l)cr“'=+i(a;fc+i)... a°'’"(xm)} = 0. 

(42) 

The hrst order equation for the Liouvillian would be given by continuing the expansion 
of Eq. (29) in lattice spacing a = 1/n, using Eq. (37) (for the Hamiltonian) and the 
explicit forms of the NESS expansion given by Eqs. (23) and (24), 


iadi/("V‘“’ + »p‘"’ =0- 


(43) 


However, we do not need this as we already know the algebraic relations of the ASOs 
dehning the steady state (Eq. ([30|) which fully determine the leading order via the 
partial differential equations ( p^ up to the boundary conditions for these equations. We 
have found these boundary conditions for the m-point correlators using the preceding 
discussion and they are given in Eqs. (40), (41) and (42), and thus we hnd. 


C^°^(T) = A]^sin(afc^ + nxk), 


k=l 


(44) 


where oa, = 0 if the operator depending on fc—th coordinate jk is or (i.e., = x, y) 

and oa; = 1 if it is cr^ (or ak = z) and A is some yet to be determined constant. 
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Let us recap: we have used the expansion of the equation for the NESS (29) in 


lattice spacing a = 1/n and via Eq. (25) found the boundary conditions for the NESS 


of an open XXZ spin chain under boundary driving given by Eq. (28). It is important 


to note that this preceding discussion on the perturbation expansion in a (via Eq. (25)), 
which we used only to hnd the boundary conditions for the partial differential equations 


(32), is not necessary provided that we also know the action of the representation of the 
ASOs (which dehne the MPA for the NESS in Eq. (|^) , i.e., we can use how the ASOs 
act on the boundary vectors |R)(|L)) to hnd the boundary conditions. This is in fact 
known for the maximally driven open XXZ spin chain [8| and some of the boundary 


conditions we found here in Eqs. (40), (41), (42) were found in the same article. 


Furthermore, in order to hnd the scaling factor A we need to employ the boundary 
conditions given by the appropriate representation of the ASOs |H]. Doing this for fairly 
large system sizes we hnd the following behaviour; if we let px denote the number of 
and Py the number of operators in the correlator then A = 0 if Pa, or py are odd, 
otherwise A is given by the following recurrence relation {p '. = Px + Py) [5T] . 




(2p)! 

u 1 = <! 


^,Px -Py-l)- A{p,Px -Py 


Px 

1)) Px 


Py=D 
Py ^ 0 


We have not been able to prove this in general, however exact calculations for up to 
u = 20 support this conjecture. The initial terms A(p, 0) and A{p + 1,0) can be 
calculated using the hrst line of the recurrence relation and, in turn, can be used to 
compute A{p +1,1), A{p + 1, 2),.... Also we take A(0, k) := 0. 

In order to calculate with higher precision than C'{n^) we will again need to 
employ the boundary conditions given by the appropriate representation of the ASOs [H]. 
These corrections will be important when we want to look at the connected correlation 
functions. For higher precision in 1/n we will likewise need a more precise scaling of the 
non-equilibrium partition function found in [^, 




6^ 



(1+<?(«-=')). 


(45) 


Interestingly, the constant a can be found through self-consistency conditions imposed 
by the algebra (30) and the boundary conditions, as will be shown later. We hnd, using 
Taylor series expansion from Eq. (18) for Eq. (30), that the Ijn corrections are 
given by the following coupled partial differential equations (omitting A = 1, setting 
fl = 4(1 — a) and switching to a more compact PDE notation). 


c+, + = (46) 

X y ((^ - 4)C'“>(5) + 2(1 - x,)C(»' + - (1 - 2x,)C<»>) + 2C+^. 

p<k,j>k 


Likewise we may hnd the higher order corrections. The PDFs determining them become 
exceedingly complicated and we omit writing them explicitly. 
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Note that these correction terms actually depend on the whether we expand the 
series in 1/n or l/(?7, — 1) (not the definition of Xk directly). The k-i\i order corrections 
are then of either order in 1/n^ or l/{n — 1)*^. Naturally, the leading orders l/n° and 
l/(n —1)° are equal, but the higher orders are not. This merely reflects the fact that the 
higher orders give increasingly precise corrections for large, but still finite, system size 
n and it is merely a matter of convention in what we will expand - they gives equivalent 
and consistent information about the correlation functions. 

We can now turn to computing the connected correlation functions. They are given 
by the standard generating function. 


~ d d 

U[ai a2 . . . J — (72 ... (7^ ) .— „ „ 


log trpoo (exp ^ Zjcr"") | 


dzi dZn 

J 

We will write out the two and three-point connected correlators explicitly. 


2,=0 


C{a 




) = (' 




C'( 


31 32 33 ' 




\ n / \ ?i / 


- cr 


)(a“V-) + 2«)(a-)(a: 


« 3 \ 


(47) 


(48) 


The expressions for higher connected correlators are quite long so we omit writing them. 


Now using Eq. (47) and the previous results from Eq. (46) (as well as the higher order 


equations obtained from Eq. (18)) and the boundary conditions given by the matrix 
representation of the ASOs pm we can find arbitrary m-point connected correlators. 
We will show how to obtain some of the lower m-point connected correlators. We 


first note that, like the correlators Eq. (44), all the connected correlators which contain 


an odd number of and operators are 0. 

First let us then start with the 2-point connected correlators; was already 

found in pm- It was shown to scale inversely with the system size, ~ 1/n. Since the 
expectation values of the transverse operators (aj) = ((tJ) = 0 we immediately see that. 




31 32 


) = a^j ) = Jsin(7rxi)sin(7ra;2), xi,2 = / , (49) 

2 77 , — 1 


where we have taken xi ^2 = to conform with previous approaches mpso]. One 

sees then that the transverse 2-point connected correlators do not decay with system size 
making them truly long-range in the sense of Ref. mi. Since the basis of decoherence 


is determined by the dissipator Eq. (28) to be in the z-direction, this can be understood 
as a purely quantum effect. 

Furthermore, when we find the I/77 correction we require also (as is given by the 
explicit representation of the ASOs pm) that the boundary conditions. 


a: 


31 J 2 7 


dxi 


I 1 1 = 

1X1 = 1,012 = 1 


(50) 
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n (cr''(xi)a>'(x2)>c 


n(cr''(xi)cr''(xi + l)a^(x2)>c 



(a''(xi)a)'(x2+i)a^(x2)>c 
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Figure 1. The two and three-point connected correlation function for A = 1 and 
e = 10. 


are satisfied. This fixes uniquely that fl = 1, or a = 3/4 in Eq. (45). Thus, 

8 

) = 2(1 - 2 / 2 ) cos(7r?/2)[7rt/i cos(7rt/i) - sin(7rt/i)] - 2yi cos(7rt/i) 


TT 


- 7r[t/i(t/i - 1) + 2 / 2(1 - 2 / 2 )] sin(7r2/i) sin(7r2/2), 


(51) 


where, as before, 0 : 1,2 = ~ min(a:i, 0 : 2 ), 2/2 = niax(a:i, 0 : 2 ). 

We hud that the mixed transverse 2-point correlator decays 

with system size as 1/n in leading order, 

= ^sin(7r(a:i - 0 : 2 )), = 0, 2 : 1,2 = ^ • (52) 

The m-point connected correlators however become quickly more complicated for higher 
orders. 
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The simplest non-zero three point connected correlation fnnction is 
scales as 1/n and is given by, 


- l)cos(7r2/3)j/2COs(7ra;i)cos(7rj/2)- 


— ^7r(xi — 1) cos^nys) [(cos(7ra:i) -f tixi sin(7rxi)) sin(7r|/2)] 

-4cos(7r?/2) [7r^|/2 cos(7ra;i) - (2-f-7r^a;i(-l 2 / 2 )) sin(7ra;i)] 

+ {8 - Stt - 27r^ [(1 - xi)xi + (1 - y2)y2 + (1 - 2/3)l/3]} cos(7ra:i) 

-h 27r^(l - 2xi) sin(7rxi) sin(7r|/2) sin^nys) , Xi < y 2 < 2 / 3 , 


(53) 


= — |27rcos(7r|/3)[-27r^a;i(-l -1- 1 / 3 ) sin(7ra;i) sin(7r|/2) 

4- cos(7rxi) cos(7r|/2) {27r^|/2(-l + ys) + (-16 + Sn - 2'Kyflj sin(7r?/2)}] 

— cos(7r|/2)(2/2 cos(7ra;i) -|- 7rxi(—1 -|- 1 / 2 ) sin(7ra;i)) 

- [(8 - Stt -4 27r^ {(-1 -f- xi)xi + (-1 + 2/2)l/2 + (-1 + 1 / 3 ) 113 )} cos(7ra;i) 

+ 2(1 + ir^(l - 2a:i))]siii(irii)sin(irj/2)siii(ir93)J, yi < x, < y, (54) 


= ~ ^i) cos(7r?/3) [nxi sin(7ra;i)) sin(7r?/2) — T^y2 cos(7ra;i) cos{ny 2 ) — (cos(7ra;i)] 

Y^l -4cos(7r|/2) [7r^|/2 cos(7rxi) - (2-4 7r^xi(-l-4 7/2)) sin(7ra;i)] 

- {8 - Svr - 27r^ [(1 4- - (1 - y 2 )y 2 + (1 - 2 / 3 )l/ 3 ]} cos(7rxi) 

-27r^(l-4 2a;i)sin(7rxi)sin(7r?/2)sin(7r|/3)i, y2<y3<Xi, (55) 


where Xk = and 1/2 = min(x 2 ,a; 3 ), y^ = max(a; 2 , 2 : 3 ). 

The simplest fonr-point fnnction is ~ c again 

long range in the sense that it doesn’t decay with system size n in the leading order, 


~ 3 jk 

(T^ cr^ cr^ ) = -- sin(7ra;i) sin(7ra;2) sin(7ra;3) sin(7ra;4), Xk = 


n — 1 


(56) 


The other correlation fnnctions stndied and plotted in Figs. and inclnde 


^ ^ and (, 




) (X 1/n. Based on onr 


resnlts, althongh we were nnable to prove so in general, we conjectnre that in leading 
order of 1/n, m-point connected correlation fnnctions containing only operators 
scale as l/n”^“^ and that all pnrely transverse m-point connected correlators, i.e., those 
containing only an even nnmber of (cr^) do not decay in the leading order, ~ (l/n)°. 
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+^)a^{x2)cf^{x2+ i))c 

0 


(cr^(xi )cr^(xi + l)cr>'(x2)cr>'(x2+ l)>c 
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Figure 2. Some of the 4-point connected correlators calculated at strong coupling 
e = 10 and A = 1. 


4 . 2 . The easy-plane regime A < 1 

For A < 1 only the cubic relation for Oz is known [39l |16] and it is given by, 

Koi'j, s)(OoOoOz -h OzOqOo) + OqOzOo + Ki(7, s){Oz, Oz} -h /t2(7, ■s)Oz = 0, (57) 

where, 

i^o{' 7 ,s) = ^ - cos ( 27 ), 

/ti( 7 , s) = 1 -|- cos ( 27 ) -f cos ( 47 ) — 4 cos ( 27 s), (58) 

^ 2 ( 7 , s) = 12 cos ( 27 s) — 2 cos ( 47 ) — 10 — 16 cos ( 27 ) sin^ ( 7 s) -|- (8 — 4 cos ( 27 s))[s]^. 
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where A = 005 ( 7 ), Q = exp(i 7 ), [x\q := — q ^)/{q — q and £[s]q = 4 icos( 7 s). 

It leads to, 

= 0, (59) 

for correlators containing only a^. This is incorrect for e.g., 2-point correlators when 
xi = X2 since = 1. This obviously means that there is a discontinuity and hence 
the assumption that we may exapand in a Taylor series, which we made when deriving 
Eq. (18), is possibly not valid. In fact computing the correlators when the difference 
between xi and X2 is of order of lattice spacing using the known discrete solution [ 8 ] 
shows that the trivial result Eq. (59) is not valid there. We cannot therefore make any 
claims near these points, but calculating these correlators for large but finite n shows 
they are simply constant functions in Xk- For finite distances among all the Xk and Xk' in 
the continuum limit (meaning that \jk — jk'\ —t 00 as n —>■ cx)) the trivial result Eq. (59) 
holds. In other words, on the ’large scale’ there is no non-smoothness. This is also the 
reason why the ’large-scale’ results also hold for A = 1 even though there (o'")((j“) = 1, 
when i = j, as well. 

Even though we could not find any similar cubic relations for 0+ nor 0_ we find 
that correlators containing these operators behave like those with only i.e., they are 
vanishing. In any case, the behavior in this regime is quite trivial as all expectation 
values of non-local operators decay to 0 exponentially fast, whereas the (ultra)-local 
expectation values are constant function throughout the chain [ 1 ]. 

For A > 1 the non-equilibrium partition function (45) decays exponentially fast 
m and we can not apply our method as all orders in l/n contribute when calculating 
the continuum limit Eq. (18). 


5. Fluctuations of spin cnrrent and magnetization operators 

Quantum mechanics is an intrinsically probabilistic theory, therefore in order to achieve 
a deeper understanding of any physical system we are studying we need to go beyond 
calculating merely the expectation values of observables we are interested in. 

We will study the probability distribution of the total magnetization in the 
maximally driven open XXZ spin chain, P{M := (t|). 

We also study the fluctuations of the instantaneous spin current flowing through 
the system once the system has reached NESS. That is, avoiding certain ambiguities 
with defining current statistics [32], we want to merely study simply the moments of 
the instantaneous spin current measured in the long time limit of the system (NESS), 
i-e., or)- 

The local spin current operator is defined via the continuity equation, = 
i[H, al] = jk — jk-i, and for the XXZ spin chain is, 

jk = 2i(cT+afc+i - (60) 

However, we notice a problem immediately, namely j™ oc jk when m is odd and 
jjj oc (1 — cr^cr^+i) when m is even which makes studying this quite trivial. 
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A possible solution to this is that we may define a multiple site current operator 


averaged over K sites, similar to the one studied in 


Ylm=k jm/K. For example, 4^^ = jk, Jlf’ = {jk + jfc+i)/2, etc. The same type of 
’space-integrated’ current was previously studied in different settings in Refs. [53l [M] • 
These operators physically correspond to the average flow of magnetization between 
site k and site K + k, i.e, 

1 + O-fc+i CtI^k) _ j{K) j(K) ^ 

K dt ~ ^ ^ 

Note that the expectation values (4^^^) = (4^^^) Ki,K 2 , which follows from 

the continuity equation in the long time limit. However, the higher moments are not 
equal for different Ki and K 2 . 

Also, these operators still have the property that with their {K + l)-th power, 
(4^^)^'*'^ oc j\^\ so we will define an extensive (up to a prefactor of l/(n — 1)) quantity 
J = where n is the system size. That is. 


r(2) 


and in 


j; 


(K) _ 


n—1 


^ = E 


Jk 


k=l 


n 


(62) 


We also define a moment generating function for the total magnetization operator, 
M = X]fc=i ^k steady state, 

(63) 


Gm(x) = (e'*") = 




6. Fluctuations of spin cnrrent and magnetization in the maximally 
bonndary driven XXZ spin chain 

Using the results from Sec. we can study the fluctuations of the spin current and 
magnetization operators, as defined in the previous section. Sec. Recall that we 
defined the averaged spin current operator as J = The ASO corresponding 

to the local spin current W ;= i(0+0_ — 0_0+) is proportional to Oq |H1 06], 

W = -2i[s]^Oo, (64) 

where, as before, A = cos( 7 ), q = exp(i 7 ), [x\q := — q~^)/{q — q~^) and 

£[s]g = 4icos(7s). The fact that W cx Oq guarantees the validity of the continuity 
equation for the magnetization,^ = i[H, al] = jk — jk-i, in the long time limit when 
the system has reached the NESS, i.e., = 0 = (jk) — (jk-i) (in other words the 

expectation values of spin current jk are equal on all sites k). 


We square (take the third power of) Eq. (62 ) and then using the definition Eq. (60), 


the properties of the Pauli matrices and the previously mentioned fact that the NESS 
expectation values (jk) = (jm) for all k,m we find (grouping equal terms together), 

{J") = (65) 


{n-iy 


{n - 2)(n - 3)(ji,2j3,4) + 


n — 1 


n—1 

E 

k=l 


WWfe+1 


+ 2(a+cT;-+2 + 
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and, 


(J^) = 


(n 


1)= 


(n - 3)(n - 4)(n - 5) (ji, 2 ^ 3 , 444 , 5 ) - 3(n - 2)(n - 3)(ji,2) 


. 5 ^ (3(a+cT;.+3 - cTfc CT++ 3 ) - (cx+a^+ia^+ 2 ^fc+ 3 ) + {a,^ 


+ “ 3)(ji,2)(cxfca^+i) - 12(n - 4)(ji,2)(cT+afc+2 + ^^^+ 2 )) 


( 66 ) 


Let us first consider the isotropic case A = 1 (g —)■ 0). The second moment (J^) contains 
sums of expectation values of three different types of operators which are non-trivial in 
the thermodynamic limit n —)■ 00 : (ji, 2 j 3 , 4 ), {^k^k +2 '^^k^k+ 2 )- Using 


Eq. (64) and the asymptotic form of the non-equilibrium partition function Eq. (45), 
it is easy to show that (ji, 2 j 3 , 4 ) oc 1 /n^ in the leading order. 

Furthermore we note that the leading order of the other terms cancels exactly. In 
fact, using, the next to leading order of which was shown in Ref. [1] to be 

(where, as before xi ,2 = 

TT r 

cos(7ra:i) 7r{(-l xi)xi + (-1 -f X 2 )x 2 } cos{'kx 2 ) 

-|- (1 — 20 : 2 ) sin(7ra;2) + sin(7ra:i){(l — 2xi) cos{'kx 2 ) — 27ra:i(—1 -f- 1 / 2 ) sin(7ra:2)} , (67) 


together with Eq. (51), we see that the next-to-leading order 1/n cancels as well. 


However expanding around X 2 = xi + 1/{n — 1) for large n we are left with a hnite 
contribution of order 1/n^. Similarly for (J^) we are left with a leading order of 1/n^ 
Finally we hnd that. 


lim (T^) = 


2 + 57r2 
8n^ 


lim (J3) = 




en° 


( 68 ) 


For a dense set of rational anisotropies, A = cos(7rm/n) with m,n G Z, it is 
known that the spin current is ballistic (reaching a constant value as n —)■ cxd) as 
the matrix representation of Oq is truncated to a hnite dimension. Thus the non- 
equilibrium partition function is determined by the largest eigenvalue of this truncated 


transfer matrix in the limit n —?• 00 HE]. Therefore in Eqs. ( [65| and ( [M| ) the terms 
(n —2)(?7, —3)(ji,2j3,4) and — (n —3)(n —4)(?7, —5)(ji,2j3,4j4,5), respectively, will dominate. 
It is clear that an analogous result holds for all the higher moments. If we let Aa denote 


the largest eigenvalue of the truncated Oq at anisotropy A then from Eq. (64) it follows. 


in = I / [s], 


(69) 


where, as always, A = cos( 7 ), q = exp(i 7 ), [x\q := (g^ — g ^)/{q — q and 
£[s]g = 4icos(7s). As s is purely imaginary [39l 06] Eq. (69) is always real. 


Finally we note that, from the dehnition Eq. (63), Eq. ( pT| and Eq. ([^, the 
moment generating function Gm{x) fh® total magnetization distribution can be 
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Figure 3. The current moments for (a) A = 1 and (b) A = 1/2 for two coupling 
parameters e = 1 and e = 10 obtained by explicit computation of the MPA for the 
NESS found in Ref. [5]. 
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Figure 4. Asymptotic scaling of the moment generating function of the total 
magnetization operator for (a) A = 1, (b) A = 1/2 obtained by explicit computation 
of Eq. (70) for the NESS using the MPA found in Ref. [5] 


given as 


Gm{x) — 


iG(xrio) 


cos(x)Oo + isin(x)0^)”|0) 


(70) 


We were unable to find closed form expression for the probability distribution of 
the total magnetization. However, we perform explicit computation of Eq. (70) using 
the MPA for NESS from Ref. [8] for large n and plot it in Fig. (|^, where we notice that 
the width of the peaks around kir seem to scale as 1/y/n. If so, this suggests that if 
we Fourier transform Eq. (70) before we take the leading order in 1/n we find that the 
probability distribution is actually a Gaussian distribution around M = 0 with variance 
OC Ij yjn. 
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Conclusions 


In this article we have elaborated on a general method to derive all the correlation 
functions of exactly solvable boundary driven quantum spin chains in the continuum 
limit. This method is based on solving partial differential equations for scaled correlators 
one obtains directly provided that a cubic algebra is satisfied by the auxiliary space 
operators forming a matrix product steady state, and incorporating also the appropriate 
boundary conditions. We utilise this method to compute various quantities for the 
exemplar case of the open boundary driven XXZ spin chain [8]. 

More specihcally, we computed all the (transverse and longitudinal) spin correlation 
functions for the boundary driven open isotropic XXX spin chain in the leading order of 
the continuum limit up to a scaling factor (for which a recursion relation was conjectured 
based on a known discrete solution [8]). We derived explicit expressions for certain (up to 
4-point) connected correlators and observed that the connected correlators of operators 
transversal to the basis of decoherence in the 2 ;-direction (i.e., tensor products of an even 
number of and cr^) exhibit long-range order — that is they do not decay with system 
size in the leading order. Similar behavior has been previously observed numerically for 
a related system im. 

Finally, we dehned two statistical quantities of interest — fluctuations of the total 
spin current (studied previously in Refs. [53l El] in different contexts and also related to 
a quantity studied in [27] and in [52]) and the magnetization operators. We computed 
the second and third moment of the spin current operator at the isotropic point A = 1 
and all the moments for A < 1 when A can be expressed via rational multiples of vr 
as A = cos(7m/m). At the isotropic point we found that for system size n the second 
moment decays oc 1/n^ and the third moment oc 1/n®. For A < 1 we find that none of 
the moments decay with the system size. 

We need to stress that it is possible that the continuum limit is not well-dehned 
when the difference between the operators in the basis for the steady state is of order of 
the lattice spacing. In this article we focused on the continuum limit of, the previously 
mentioned, open XXZ spin chain for which a discrete solution is known [8] and this 
allowed us to check the validity of the continuum limit. 

Even though our treatment for the continuum limit is exact, it would be interesting 
to see whether one can End a solution for an open quantum non-equilibrium steady 
state directly for a continuum system (and thus avoiding knowing an exact solution for 
any hnite, discrete system), either on the level of a quantum field theory or at least 
perturbatively in lattice spacing in the leading order (in the sense of Eq. (25)). 
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